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Discretization of the viscous terms in current finite-volume unstructured-grid schemes 
are compared using node-centered and cell-centered approaches in two dimensions. Ac- 
curacy and complexity are studied for four nominally second-order accurate schemes: a 
node-centered scheme and three cell-centered schemes - a node-averaging scheme and 
two schemes with nearest-neighbor and adaptive compact stencils for least-square face 
gradient reconstruction. The grids considered range from structured (regular) grids to 
irregular grids composed of arbitrary mixtures of triangles and quadrilaterals, includ- 
ing random perturbations of the grid points to bring out the worst possible behavior of 
the solution. Two classes of tests are considered. The first class of tests involves smooth 
manufactured solutions on both isotropic and highly anisotropic grids with discontinuous 
metrics, typical of those encountered in grid adaptation. The second class concerns solu- 
tions and grids varying strongly anisotropically over a curved body, typical of those en- 
countered in high-Reynolds number turbulent flow simulations. Tests from the first class 
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indicate the face least-square methods, the node-averaging method without clipping, and 
the node-centered method demonstrate second-order convergence of discretization er- 
rors with very similar accuracies per degree of freedom. The tests of the second class are 
more discriminating. The node-centered scheme is always second order with an accuracy 
and complexity in linearization comparable to the best of the cell-centered schemes. In 
comparison, the cell-centered node-averaging schemes may degenerate on mixed grids, 
have a higher complexity in linearization, and can fail to converge to the exact solu- 
tion when clipping of the node-averaged values is used. The cell-centered schemes using 
least-square face gradient reconstruction have more compact stencils with a complexity 
similar to that of the node-centered scheme. For simulations on highly anisotropic curved 
grids, the least-square methods have to be amended either by introducing a local map- 
ping based on a distance function commonly available in practical schemes or modifying 
the scheme stencil to reflect the direction of strong coupling. The major conclusion is 
that accuracies of the node centered and the best cell-centered schemes are comparable 
at equivalent number of degrees of freedom. 


I. Introduction 

Both node-centered and cell-centered finite-volume discretizations are widely used for com- 
plex three-dimensional turbulent simulations in aerospace applications. The relative advantages 
of the two approaches have been extensively studied in the search for methods that are accurate, 
efficient, and robust over the broadest possible range of grid and solution parameters. The topic 
was discussed in a panel session at the 2007 AIAA Computational Fluid Dynamics conference, but 
a consensus did not emerge. One of the difficulties in assessing the two approaches is that compar- 
ative calculations were not completed in a controlled environment, i.e., computations were made 
with different codes and different degrees of freedom and the exact solutions were not known. 

In this paper, we provide a controlled environment for comparing a subset of the discretization 
elements needed in turbulent simulations, namely that of the viscous discretization. In particular, 
we consider only Poisson’s equation as a model of viscous discretization. We use the method of 
manufactured solution, so that the exact solution is known, and conduct theoretical and computa- 
tional studies of both accuracy and efficiency for a range of grids. 

The two-dimensional (2D) grids considered range from structured (regular) grids to irregular 
grids composed of arbitrary mixtures of triangles and quadrilaterals. Highly irregular grids are 
deliberately constructed through random perturbations of structured grids to bring out the worst 
possible behavior of the solution. Two classes of tests are considered. The first class of tests 
involves smooth manufactured solutions on both isotropic and highly anisotropic grids with dis- 
continuous metrics, typical of those encountered in grid adaptation. The second class of tests 
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concerns solutions and grids varying strongly anisotropically over a curved body, typical of those 
encountered in high-Reynolds number turbulent flow simulations. 

There are four main schemes considered, a node-centered (NC) scheme and three cell-centered 
(CC) schemes. The cell centered schemes include a node- averaging (CC-NA) scheme and two 
least-square face gradient reconstruction schemes differing in their stencils - a nearest-neighbor 
(CC-NN) stencil and an adaptive compact (CC-CS) stencil. The effect of clipping is studied for the 
CC-NA scheme. In its current version, the CC-CS scheme is applied only on triangular grids. For 
the second class of tests, an approximately mapped (AM) least-square approach is introduced to 
accommodate curved high-aspect-ratio grids. The mapping uses the distance function commonly 
available in practical codes and can be used with any scheme. Each of the schemes considered 
is nominally second-order accurate. The properties to be compared are computational complexity 
(operation count) and discretization accuracy at equivalent degrees of freedom. 

II. Grid terminology 

This paper studies finite-volume discretization (FVD) schemes for viscous fluxes on grids that 
are loosely defined as irregular. There is no commonly accepted definition for irregular grids, so, 
for clarity, we specify the grid terminology used in this paper. 

A grid is classified as periodic if it has (1) a periodic node connectivity pattern (i.e., the num- 
ber of edges per node changes periodically) and (2) a periodic cell distribution (i.e., the grid is 
composed of periodically repeated combinations of cells). Thus, periodic grids can be analyzed 
by Fourier analysis. Grids that are derived from periodic grids by a smooth mapping are called 
regular grids. Regular grids include, but are not limited to, grids derived from Cartesian ones - 
triangular grids obtained by diagonal splitting with a periodic pattern, smoothly stretched grids, 
skewed grids, smooth curvilinear grids, etc. Grids that cannot be smoothly mapped to a periodic 
grid are called irregular grids. Grids with varying local topology are called unstructured , e.g., 
grids with the number of edges changing from node to node with no pattern. 

The regular and irregular grids considered in this paper are derived from an underlying (possi- 
bly mapped) Cartesian grid with meshsizes h x and h y and the aspect ratio A — both meshsizes 
of the underlying grid are assumed to be small, h y < 1 ,h x < 1. Irregularities are introduced 
locally and do not affect grid topology and metrics outside of a few neighboring cells. A local 
grid perturbation is called random if it is independent of local perturbations introduced beyond 
some immediate neighborhood. For computational grids generated for the reported studies, grid 
irregularities are introduced in two ways (both local and random): (1) the quadrilateral cells of 
the underlying grid are randomly split (or not split) into triangles; (2) the grid nodes are perturbed 
from their original positions by random shifts; the shifts are fractions of a local meshsize. 
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Six grid types are considered: (I) regular quadrilateral (i.e., mapped Cartesian) grids; (II) 
regular triangular grids derived from the regular quadrilateral grids by the same diagonal split- 
ting of each quadrangle; (III) random triangular grids , in which regular quadrangles are split by 
randomly chosen diagonals, each diagonal orientation occurring with probability of half; (IV) per- 
turbed triangular grids , which are random triangular grids, with grid nodes perturbed from their 
initial positions by random shifts; (V) perturbed quadrilateral grids, which are quadrilateral grids 
with randomly perturbed grid nodes; and (VI) perturbed mixed-element grids, in which perturbed 
quadrangles are randomly split or not split by randomly chosen diagonals. Grids of types (III), 
(IV), and (VI) are irregular because there is no periodic connectivity pattern. Grids of types (IV)- 
(VI) are irregular because there is no periodic cell distribution. The representative grids are shown 
in Figure 1 . 

Our main interest is the accuracy of FVD schemes on general irregular (mostly unstructured) 
grids with a minimum set of constraints. In particular, we do not require any grid smoothness, 
neither on individual grids nor in the limit of grid refinement. The only major requirement for a 
sequence of refined grids is to satisfy the consistent refinement property. The property requires 
the maximum distance across the grid cells to decrease consistently with increase of the total 
number of grid points, N. In particular, the maximum distance should tend to zero as iW 1 / 2 in 2D 
computations. For three-dimensional (3D) unstructured grids, the consistent refinement property 
is studied in [1]. On 2D grids, the effective meshsize, h e , is computed as the L 1 norm of the square 
root of the control volumes. 

The locations of discrete solutions are called data points. For consistency with the 3D ter- 
minology, the 2D cell boundaries are called faces, and the term “edge” refers to a line, possibly 
virtual, connecting the neighboring data points. Each face is characterized by two vectors: (1) the 
edge vector, which connects the data points of the cells sharing the face and (2) the directed-area 
vector, which is normal to the face and has the amplitude equal to the face area. For each cell/face 
combination, the vectors are directed outward. Local grid skewness is measured at a face as the 
skew angle between the edge vector and the directed-area vector; small skew angles correspond to 
low skewness. On grids with convex cells, the maximum skewness is bounded by 90°; on more 
general grids with non-convex cells, the skew angle can approach 180°. 

For most of the computational tests, the random node perturbation in each dimension is defined 
as \rh, where r G [—1, 1] is a random number and h is the local meshsize along the given dimen- 
sion. With these perturbations, triangular cells in the rectangular geometry are allowed to collapse, 
i.e., a cell may become a zero-volume cell, albeit with a probability approaching zero. The random 
perturbations are introduced independently on all grids in grid refinement implying that grids of 
types (IV)-(VI) are grids with discontinuous metrics, e.g., ratios of neighboring cell volumes and 
face areas are random on all grids and do not approach unity in the limit of grid refinement. 
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(c) Type (III): Random triangular grid 


(d) Type (IV): Perturbed triangular grid 




(e) Type (V): Perturbed quadrilateral grid (f) Type (VI): Perturbed mixed grid 

Figure 1. Typical regular and irregular high-aspect-ratio grids with vertical scale expanded for better visual- 
ization. 
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III. Finite- volume discretization schemes 


The considered model problem is the Poisson equation 


A U = f, 


( 1 ) 


subject to Dirichlet boundary conditions; function / is a forcing function independent of the solu- 
tion. The 2D primal meshes generated for this study are composed of triangular and quadrilateral 
cells. The FVD schemes are derived from the integral form of a conservation law 



( 2 ) 


where VC/ is the solution gradient, 0 is a control volume with boundary T, and n is the outward 
unit normal vector. The general FVD approach requires partitioning the domain into a set of 
non-overlapping control volumes and numerically implementing equation (2) over each control 
volume. 

Cell-centered and node-centered discretizations are con- 
sidered. Cell-centered discretizations assume solutions are 
defined at the centers of the primal grid cells with the pri- 
mal cells serving as the control volumes. The cell center 
coordinates are typically defined as the averages of the co- 
ordinates of the cell’s vertexes. Node-centered discretiza- 
tions assume solutions are defined at the primal mesh nodes. 

For node-centered schemes, control volumes are constructed 
around the mesh nodes by the median-dual partition: the cen- 
ters of primal cells are connected with the midpoints of the 
surrounding faces. These non-overlapping control volumes 

cover the entire computational domain and compose a mesh Figure 2 ‘ Control-volume partitioning 

for finite-volume discretizations. Num- 

that is dual to the primal mesh. Both cell-centered and node- , . iri ^ . r , 

r bers 0—12 and letters A — L denote grid 

centered control-volume partitions are illustrated in Figure 2. noc | es and primal cell centers, respec- 
tively. The control volume for a node- 
A. Cell-centered FVD schemes centered discretization around the grid 

node 0 is shaded. The control volume for 
In cell-centered discretizations, the conservation law (2) is a cell-centered discretization around the 
enforced on control volumes that are primary cells. The flux cell center A is hashed, 
at a face is computed as the inner product of the solution 
gradient at the face and the directed-area vector. With refer- 
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ence to Figure 2, the gradient, VC/o 4 , at the face linking nodes 0 and 4 is found by satisfying two 
relations: 

1. (VC / 0 4 • e AB C ) — rf-^TT’ where e^| e = is the edge [ AB ] unit vector; Ua and Ub 

are the solutions defined at the cell centers; and r c B are the cell-center coordinate vectors. 

2. (V!7o 4 • eQ 4 Ce ) = I Vf/^ ace |, where e{^ ce = , r ^~ r ° . is the face unit vector and I \/U^ ace \ is the 
solution derivative computed along the face [0,4]; r" is a coordinate vector of the node i. 

The cell-centered FVD schemes considered in this paper differ only in computing | Vf/^ ace |. 

1. Cell-centered node-averaging FVD schemes 

In the CC-NA schemes, the solution derivative along the face is computed as the divided difference 
between the solution values reconstructed at the nodes from the surrounding cell centers. With 
respect to Figure 2, the solution at the node 0 is reconstructed by averaging solutions defined 
at the cell centers A , B, and C. The solution reconstruction proposed in [2, 3] and used in [4] 
is an averaging procedure that is based on a constrained optimization to satisfy some Laplacian 
properties. The scheme is second-order accurate and stable when the coefficients of the introduced 
pseudo-Laplacian operator are close to 1. It has been shown 5 that this averaging procedure is 
equivalent to an unweighted least-square linear fit. The gradient is resolved from the derivative 
along the face and the derivative along the edge connecting the cell centers across the face. For the 
face [0,4], the face derivative is computed using solutions reconstructed by node averaging at the 
nodes 0 and 4; the edge derivative is computed using solutions at cell centers A and B. 

On highly stretched and deformed grids, some coefficients of the pseudo-Laplacian may be- 
come negative or larger than 2, which has a detrimental effect on stability and robustness. 6,7 
Holmes and Connell 2 proposed to enforce stability by clipping the coefficients between 0 and 

2. The CC-NA schemes with clipping represent a current standard in practical computational fluid 
dynamics (CFD) for applications involving cell-centered finite volume formulations. 8 As shown 
further in the paper, clipping seriously degrades the solution accuracy. 

2. Cell-centered face-least-square scheme 

An alternative cell-centered scheme relies on a face-based least-square method. First, an auxiliary 
face gradient is reconstructed within a face using a least-square procedure. The derivative along 
the face, |VC^ ace |, is computed by inner product of the auxiliary gradient with the face tangent 
(unit) vector. 

The two approaches to determine stencils for the least-square linear fit at a face are described 
as follows. The nearest-neighbor CC-NN six-point stencil includes the two prime cells sharing the 
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(a) Face least-square stencil with CC-NN 
scheme 


(b) Face least-square stencil with CC-CS 
scheme 
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(c) Cell Laplacian stencil with CC-NN 
scheme 
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(d) Cell Laplacian stencil with CC-CS 
scheme 


Figure 3. Stencils on high-aspect-ratio grids of type (III). Figures are vertically expanded for better visualiza- 
tion. 


face and their face neighbors that share one of the face nodes. In Figure 3(a), the CC-NN stencil 
for the highlighted face is denoted by circles. 

The adaptive compact stencil (CC-CS) is important for discretizations on high-aspect-ratio 
grids of types (II) and (III) to represent correctly the direction of the strong coupling. It is con- 
structed by choosing between two stencils for face least-square gradient reconstruction: a six-point 
stencil and a four-point stencil. In general, the four-point stencil is intended for “long” faces of 
high-aspect-ratio grids in local geometries resembling a layer formation. 

Technically, at each face, the CC-CS method first attempts to construct a six-point stencil by 
combining two prime cells and four auxiliary cells; each auxiliary cell is associated with a prime 
cell and a face node. The method chooses the auxiliary cell that (1) shares the face node, (2) is 
located on the opposite side of the face from the associated prime cell center, (3) is not already 
in the stencil, and (4) has the shortest distance to the center of the associated prime cell. The 
six-point stencil for the highlighted diagonal face is denoted by the union of circles and crosses 
in Figure 3(b). In the process of construction, each prime cell identifies the closest associated 
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auxiliary cell, which is called its mirror image. The four-point stencil is defined by prime cells and 
their mirror images and shown as circles in Figure 3(b). 

The CC-CS method selects the four-point stencil, if either the six-point stencil cannot be formed 
following the rules (l)-(4), which may happen next to the boundaries or in curved geometries, or 
the four-point stencil represents an ideal pair-wise construction. The four-point pair-wise con- 
struction is considered ideal, if the data points within the pairs of prime-cell/mirror-image are very 
close (closer than a predefined distance threshold) and have a low skewing (the angle between the 
vector connecting the prime cell center with its mirror image and the face directed area vector is 
smaller than a predefined skew threshold). For computations on high-aspect-ratio grids, the dis- 
tance threshold has been chosen as three times the smaller mesh size of the background grid and 
the skew threshold has been set to sin _1 (0.1). The four-point stencil in Figure 3(b) is considered 
ideal. 

Figures 3(c) and 3(d) compare CC-NN and CC-CS FVD stencils on the shaded cell. CC-CS 
FVD scheme uses four-point stencils for diagonal and horizontal faces and a six-point stencil for 
the vertical face. The CC-CS stencil is more compact than the CC-NN stencil and provides a three- 
point vertical structure centered at the shaded cell center that better reflects the grid anisotropy 
direction. 

Remark. It is known that on high-aspect-ratio curved grids, unweighted least- square methods 
have difficulties with reconstructing accurate gradients within a cell. 9 " 11 Inverse-distance weight- 
ing has been shown to improve gradient accuracy. For face-centered least-square reconstruction, 
the usual weightings, with distances measured from the face center, do not improve gradient accu- 
racy because all points involved in least-square stencils are typically at comparable distances from 
the face center. A modified weighting, which is based on (minimal) distances from the two cell cen- 
ters across the face, with an extended stencil (the stencil that is used in CC-NA scheme) improves 
gradient accuracy on high-aspect-ratio curved grids derived by an advanced layer method. The 
weighting effectively reduces the extended stencil to the four-point stencil of the CC-CS scheme. 
However, the method led to unstable formulations on general irregular grids and was not pursued 
further. 

B. Node-centered FVD scheme 

The second-order accurate node-centered FVD scheme illustrated by Figure 4 represents a standard 
CFD approach to node-centered viscous discretizations. The scheme approximates the integral flux 
through the dual faces adjacent to the edge [P 0 , P 4 ] as 

f VU-ndT = VU R -n R + VU L -n L . (3) 

ABC 
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The gradient is reconstructed separately at each dual face as follows. 

For the triangular element contribution, the gradient is 
determined from a Green-Gauss evaluation at the primal-grid 
element, 

VU L = W014. (4) 

p 2 

The gradient overbar denotes a gradient evaluated by the 9 
Green-Gauss formula on the primal cell identified by the 
point subscripts. With fully-triangular elements, the formu- 
lation is equivalent to a Galerkin finite element scheme with 
a linear basis function. 6 - 12 Analysis in Appendix A shows 
that on triangular grids of types (II) and (III) in rectangu- 
lar geometries, the formulation recovers the five-point Lapla- 
cian stencil of the type (I) grids, independent of aspect ratio. 

For the quadrilateral element contribution, the gradient is evaluated as 


• P ” 


Figure 4. Illustration of gradient re- 
construction for viscous terms on mixed 
grids with median-dual partition. 


VU R = W0234 + NVEri-Vtf, 


'0234 


604] eo 4 , 


(5) 


where 


I vtn = 


Ua — Un 


n - If 


(6) 


is the edge gradient computed at the median Pb of the edge [P 0 , P 4 ], U t and r? are the solution and 
the coordinate vector of the node Pj, and 


©04 — 


1, - If 


If - In 


(7) 


is the unit vector aligned with the edge [P 0) A]- Note that for grids with dual faces perpendicular 
to the edges, the edge-gradient, \'VU E \, is the only contributor. This approach to the gradient 
reconstruction is used to decrease the scheme susceptibility to odd-even decoupling. 13 - 14 It has 
been shown 115 that the scheme possesses second-order accuracy for viscous fluxes on general 
isotropic mixed-element grids. 


IV. Complexity of discretization stencils 

The size of the stencil for the viscous discretization is examined for 2D and 3D cell-centered 
and node-centered FVD schemes. Estimates are made for Cartesian meshes split into triangular and 
tetrahedral elements and neglecting any boundary effects. Estimates are compared to a numerical 
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calculation on an actual 3D turbulent viscous grid that includes boundary effects. 

In two dimensions, two splittings of the Cartesian grid are considered. The first splits each 
quadrilateral cell with a diagonal oriented in the same direction. The second splits the cells with 
diagonals of face-adjacent quadrilaterals oriented in the opposite direction. The second splitting 
is slightly more analogous to the 3D splitting. In three dimensions, half of the grid nodes have 
18 incident edges (32 incident tetrahedra) and half have 6 incident edges (8 incident tetrahedra). 
Each of the tetrahedra interior to an originally-hexahedral cell is defined by four nodes, each with 
18 incident edges. Each of the four surrounding tetrahedra within an originally-hexahedral cell is 
defined by three nodes with 18 incident edges and 1 node with 6 incident edges. For reference, 
Table 1 shows the average and maximum number of edges, n e d ge , connecting to a grid node. The 
average number of connecting edges sets the stencil size for the node-centered scheme as s = 
fi e dge + 1. The number of connecting edges is also an important factor in the CC-NA schemes 
because the gradients in the control volumes faces are computed from cell-centered data averaged 
to the grid nodes. 8 


Dimension 

n e dge (Average) 

nedge (Maximum) 

2D 

6 

8 

3D 

12 

18 


Table 1. Edges connecting to a grid node in the split Cartesian grids. 


Node-centered 

CC-NA 

CC-NN 

7 

13/16 

10/9 


Table 2. Size of the viscous stencil for 2D discretizations with triangular elements. 



Node-centered 

CC-NA 

CC-NN 

Estimate 

13 

79 

15 

Numerical 

14 

69 

15 


Table 3. Average size of the viscous stencil for 3D discretizations with tetrahedral elements. 

Table 2 shows stencil-size estimates for 2D triangular grids. There is a slight difference in 
the estimates from the two splittings (entries separated by slashes in the table), depending on the 
diagonalization pattern. The CC-NA stencil is the largest. Table 3 shows stencil-size estimates for 
3D tetrahedral grids. The CC-NA stencil is again the largest. The CC-NN stencil is only slightly 
larger than the stencil of the node-centered discretization, in both estimation and computation. The 
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complexity of the CC-CS stencil is even smaller. 

Because the computational stencil of the CC-NA scheme is very large, direct relaxation of the 
full linearization operator is prohibitively expensive in terms of CPU time and memory resources. 
Defect-correction iterations are widely used as a practical solution method. A common driver used 
for defect-correction iterations is an edge-based thin- shear-layer (TSL) FVD scheme. Although 
the TSL schemes are not accurate on general grids (accuracy of TSL FVD schemes is zeroth order 
on general grids), the defect-correction method is efficient on grids with low skewing. It has been 
shown 16 that for high grid skewing, unavoidable in highly anisotropic triangular/tetrahedral grids 
used for high Reynolds number simulations, efficiency and stability of defect-correction methods 
with TSL drivers degrade leading to explicit stability (CFL) limitations. 

V. Analysis methods 

A. Method of manufactured solution 

Accuracy of FVD schemes is analyzed for known exact or manufactured solutions. The forcing 
function and boundary values are found by substituting this solution into the Poisson equation with 
Dirichlet boundary conditions. The discrete forcing function is defined at the data points. 

1. Discretization error 

The main accuracy measure is the discretization error , E d , which is defined as the difference 
between the exact discrete solution, U h , of the discretized equations (2) and the exact continuous 
solution, U , to the differential equation (1) 

E d — U — U h ; (8) 


U is sampled at data points. 

2. Truncation error 

Another accuracy measure commonly used in computations is truncation error. Truncation error, 
E t , characterizes the accuracy of approximating the differential equation (1). For finite differences, 
it is defined as the residual obtained after substituting the exact solution U into the discretized 
differential equations. 17 For FVD schemes, the traditional truncation error is usually defined from 
the time-dependent standpoint. 18,19 In the steady-state limit, it is defined (e.g., in [20]) as the 
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residual computed after substituting U into the normalized discrete equations (2), 


E t = 


1 

W\ 



(VC/ • n) dT , 


n 


r 


where |0| is the measure of the control volume, 


(9) 


\n\ = If <jn, (io) 

n 

f h is an approximation of the forcing function / on O, and the integrals are computed according to 
some quadrature formulas. Note that convergence of truncation errors is expected to show the order 
property only on regular grids; on irregular grids, it has been long known that the design-order 
discretization-error convergence can be achieved even when truncation errors exhibit a lower-order 
convergence or, in some cases, do not converge at all. 20-22 


3. Accuracy of gradient reconstruction 

Yet another important accuracy measure is the accuracy of gradient approximation at a control- 
volume face. For second-order convergence of discretization errors, the gradient is usually required 
to be approximated with at least first order. For each face, accuracy of the gradient is evaluated 
by comparing the reconstructed gradient, V r , with the exact gradient, V exact, computed at the face 
center. The accuracy of gradient reconstruction is measured as the relative gradient error: 


-Erel — 


lie'll’ 


( 11 ) 


where functions e and G define face-wise amplitudes of the gradient error and the exact gradient, 
respectively, 

6 = \ V r U h - V exact U | , and G = |V exact (7|; (12) 

U and U h are a differentiable manufactured solution and its discrete representation (usually injec- 
tion) on a given grid, respectively; || - 1| is a norm of interest computed over the entire computational 
domain. 
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VI. Isotropic irregular grids 


A. Grid refinement 

A sequence of consistently refined grids of type (IV) is generated on the unit square [0,1] x [0,1]. 
Irregularities are introduced at each grid independently, so the grid metrics remain discontinuous 
on all the grids. The ratio of areas of neighboring faces can be as large as 3\/2; because a control 
volume can be arbitrary small, the ratio of the neighboring volumes can be arbitrary high. Isotropic 
grids randomly generated for this study have 0.01% of cell volumes smaller than -k-Tj ; the total 
number of grids nodes is ( N + l) 2 ; about 0.006% of interior faces have skew angles larger than 
70°. 

B. Gradient reconstruction accuracy 
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Figure 5. Accuracy of gradient reconstruction on isotropic irregular grids 

The accuracy of gradient reconstruction for isotropic irregular grids is first order for all meth- 
ods, 23 which is sufficient for second-order discretization accuracy. As an example, the gradient 
reconstruction tests are performed for the manufactured solution U = sin ( ttx + 2i ly). Figure 5 
shows convergence of the norms of relative gradient errors computed for cell-centered schemes 
on a sequence of refined grids of type (IV). All methods provide first-order gradient approxima- 
tions and very similar relative errors. 
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C. Convergence of truncation and discretization error 


The numerical tests evaluating convergence of truncation and discretization errors are performed 
with Dirichlet boundary conditions specified from the manufactured solution U = sin (nx + 2ny). 
For cell-centered formulations, the solution is specified on all cells linked to the boundary. Figure 6 
shows convergence of the L\ norms of truncation and discretization errors for the node-centered 
and two cell-centered formulations on grids of type (IV). As predicted in [1, 15], truncation errors 
do not converge on irregular grids in any norm. Discretization errors converge with second order 
for all formulations considered. The discretization errors of the cell-centered and node-centered 
FVD schemes are almost over-plotted, indicating a similar accuracy per degree of freedom, which 
implies better accuracy for the cell-centered FVD schemes on the same grids, as expected. 



10 2 10 1 
Effective Mesh Size 



(a) Truncation Errors (b) Discretization Errors 

Figure 6. Convergence of L \ -norms of truncation and discretization errors on random triangular grids of type 
(IV). 

D. Effects of clipping 

The tests reported in this section are performed for the CC-NA schemes and demonstrate detrimen- 
tal effects of clipping on accuracy of gradient approximation and on the discretization accuracy. 
The accuracy is evaluated for the manufactured solution U = sin(27n/). Dirichlet boundary condi- 
tions over- specified from the manufactured solution are imposed at all control volumes connecting 
to the boundary. Considered irregular grids of type (IV) are derived from underlying isotropic (unit 
aspect ratio) Cartesian grids covering the unit square. Figure 7(a) shows an example of an isotropic 
random triangular grid of type (IV) with 17 2 nodes; nodes where clipping occurs in reconstruction 
are circled; about 7% of the interior nodes are clipped. 

It has been demonstrated in [16] that the face gradients computed by the CC-NA scheme with 
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(a) Random triangular grid with 17 2 nodes. Clipped nodes are 
circled. 


E 

o 


Q 



effective meshsize, h 

e 


(b) Discretization errors 

Figure 7. Accuracy of CC-NA schemes on isotropic irregular triangular grids of type (IV) 


clipping do not approximate the exact gradients on grids of type (IV). The normal and tangential 
components of the computational gradients have been evaluated within interior faces and compared 
with the exact-gradient components at the face center. The maximum norms of the deviations be- 
tween the computational and the exact- gradient components have not converged in grid refinement. 
The CC-NA scheme without clipping provides a first-order accurate gradient approximation. Fig- 
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ure 7(b) exhibits convergence of the L \ norms of discretization errors. The no-clipping CC-NA 
scheme demonstrates second-order convergence; convergence of the CC-NA scheme with clip- 
ping appears as second order on the coarse grids, but then degrades to zeroth order. Although not 
shown, the L ^ norms of discretization errors converge with the same orders as the corresponding 
Li norms. 


VII. Anisotropic irregular grids 


0.6 r 



-0.-1 1 1 1 1 1 1 1 1 1 1 1 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


Figure 8. Random triangular stretched grid with 17 x 65 nodes 

In this section, we study FVD schemes on irregular stretched grids generated on rectangular 
domains. Figure 8 shows an example grid with the maximal aspect ratio A = 1000. A sequence 
of consistently refined stretched grids is generated on the rectangle (x, y ) € [0, 1] x [0, 0.5] in the 
following 3 steps. 

1. A background regular rectangular grid with ( N x + 1) x (N y + 1) nodes and the horizontal 
mesh spacing h x = jj- is stretched toward the horizontal line y = 0.25. The ^-coordinates 
of the horizontal grid lines in the top half of the domain are defined as 

- • ( N y iA N 

r/iVy = 0.25; y y = yj-\ + h y (3 J v 2 / , j = — ^ + 2, . . . , N y , N y + 1. (13) 

2 ' x Z 

Here h y = is the minimal mesh spacing between the vertical lines; A = 1000 is a fixed 
maximal aspect ratio; /3 is a stretching factor, which is found from the condition y,y : , ; +i = 1. 
The stretching in the bottom half of the domain is defined analogously. 
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2. Irregularities are introduced by random shifts of interior nodes in the vertical and horizontal 
directions. The vertical shift is defined as A yj = ^pmin(/i^ -1 , h J y ), where p is a random 
number between —1 and 1, and h : y ~ 1 and h J y are vertical mesh spacings on the background 
stretched mesh around the grid node. The horizontal shift is introduced analogously, A Xi = 
jQph x . With this random node perturbations, all perturbed quadrilateral cells are convex. 

3. Each perturbed quadrilateral is randomly triangulated with one of the two diagonal choices; 
each choice occurs with a probability of one half. 

A recent study 23 assessed accuracy of gradient approximation on various irregular grids with 
high aspect ratio A = ^ 1. The study indicates that for rectangular geometries and functions 

predominantly varying in the direction of small mesh spacing (//-direction), gradient reconstruction 
is accurate. For manufactured solutions significantly varying in the direction of larger mesh spacing 
(x-dircction), the face-gradient reconstruction may produce extremely large ()(Ah :r ) relative errors 
affecting the accuracy of the //-directional gradient component. Figure 9 shows examples of first- 
order accurate gradient approximations that exhibit large relative errors on high-aspect-ratio grids 
of type (III). 



Figure 9. Large relative errors in approximation of face gradients are observed for the manufactured solution 
U = sii i (tt x + 2 Try) on anisotropic grids of type (III) with A = 10 6 . 

A summary of the previous results 23 for grids of types (I)-(VI) (supplemented by the results 
for the CC-CS scheme) is presented in Table 4. All considered gradient reconstruction methods 
may generate large relative errors on irregular grids of types (IV)-(VI) with perturbed nodes. The 
CC-NA and CC-NN methods may also have large relative errors on grids of types (II) and (III); 
the NC method using the Green-Gauss approach and the CC-CS method always provide accurate 
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gradients on these grids. 


Table 4. Relative error of gradient reconstruction 


Grids 

(I) 

(II) 

(HI) 

(IV) -(VI) 

NC 

o(K 

oK 

o(h x ) 

0(Ah x ) 

CC-NA 

OK 

O(Ahl) 

0(Ah x ) 

0(Ah x ) 

CC-NN 

O(hl) 

O(Ahl) 

0(Ah x ) 

0(Ah x ) 

CC-CS 

0{hl) 

OK 

0(h x ) 

0(Ah x ) 



Figure 10. Convergence of discretization errors for solution U = sin (tt.x + 2tt y) on stretched grids of type (IV) 

A poor gradient reconstruction accuracy, however, does not necessarily imply large discretiza- 
tion error. Mavriplis 9 reported (second-order) accurate node-centered solutions even on grids with 
large gradient reconstruction errors. Here, we observe similar results for cell-centered and node- 
centered formulations. 

Convergence of the L\ norms of discretization errors for the manufactured solution U = 
sin ( 7 tx + 2rcy) on a sequence of consistently refined stretched grids of type (IV) (see Figure 8) 
is shown in Figure 10. All tests have been performed stochastically, i.e., several grids (ten) with 
different irregularities patterns have been independently generated on each scale (same number of 
nodes). The plot symbols indicate the mean errors and bars indicate the maximum and minimum 
errors observed on each scale. The effective meshsize is practically the same for all cell-centered 
schemes at a given scale, but, for visualization purposes, plots of the CC-NA and CC-CS schemes 


19 of 29 


are shifted to the right and the left, respectively, of the CC-NN scheme. 

All discretization errors are relatively small and converge with second order. The error of the 
CC-CS scheme is the smallest. The variation of discretization errors at a given scale is largest with 
the CC-NN scheme. The NC scheme is remarkably insensitive to grid irregularities. This property 
is probably explained by smaller variations (less significant irregularity) of control volumes ob- 
served for the node-centered formulation on the considered grids; the ratio of neighboring control 
volumes is always bounded and typically close to unity. As mentioned above, in cell-centered for- 
mulations on grids of types (IV)-(VI), the ratio of neighboring control volumes can be very high. 
Note that for cell-centered formulations on grids of types (I)-(III), the ratio of the neighboring con- 
trol volumes is also close to unity; thus, error variations on such grids are expected to be smaller. 
Error variations for all schemes are getting smaller on finer grids. Overall, the levels of discretiza- 
tion errors obtained with the node-centered and cell-centered schemes are similar at comparable 
degrees of freedom. 

VIII. Grids with curvature and high aspect ratio 

In this section, we discuss accuracy of FVD schemes on grids with large deformations induced 
by a combination of curvature and high aspect ratio. The grid nodes are generated from a cylin- 
drical mapping where (r, 6) denote polar coordinates with spacings of h r and h g , respectively; the 
innermost radius is r = R. The grid aspect ratio is defined as the ratio of mesh sizes in the cir- 
cumferential and the radial directions, A = The mesh deformation is characterized by the 
parameter T: 

r = fl(l-cos(h 0 )) _ Rh | = A h e _ 
h r 2 h r 2 

The following assumptions are made about the range of parameters: fi ^ 1, d > 1, and 
T h r <C 1, which implies that both h r and h g are small. For a given value of A, the parameter T may 
vary: T > 1 corresponds to meshes with large curvature-induced deformation; T C 1 indicates 
meshes that are locally (almost) Cartesian. In a mesh refinement that keeps A fixed, T = 0(Ah g ) 
asymptotes to zero. This property implies that on fine enough grids with fixed curvature and aspect 
ratio the discretization error convergence is expected to be the same as on similar grids generated 
on rectangular domains with no curvature. 

We are focusing on convergence of discretization errors on high-T grids with large curvature- 
induced deformations, following a previous study 23 focused on gradient accuracy. Considered 
manufactured solutions predominantly vary in the radial direction of small mesh spacing. 

Four types of 2D grids are considered for the cylindrical geometry: (I) regular quadrilateral 
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grids; (II) regular triangular grids derived from the regular quadrilateral grids by the same diag- 
onal splitting of each quadrangle; (III) random triangular grids , in which regular quadrangles are 
split by randomly chosen diagonals, each diagonal orientation occurring with probability of half; 
(VI) mixed grids , in which regular quadrangles are randomly split or not split by randomly chosen 
diagonals; about half of quadrangles are split. In distinction from grids of type (VI) introduced 
in the previous sections, random node perturbation is not applied to high-r cylindrical grids be- 
cause even small perturbations in the circumferential direction may lead to non-physical control 
volumes. Representative stretched grids of types (III) and (VI) are shown in Figures 11. 




(a) Grid of type (III) 


(b) Grid of type (VI) 


Figure 11. Representative 9 x 33 stretched high-F grids 


A. Accuracy of gradient approximation 

Previous observations and analysis 9-11 confirmed that the unweighted-least- square gradient ap- 
proximation is zeroth order accurate on deformed grids with high F. To improve the accuracy 
of gradient approximation, a least-square minimization in a mapped domain is proposed. Two 
mapped least-square methods are considered: an exact mapping (EM) method that uses the polar 
coordinates directly and a more general approximate mapping (AM) method that is based on the 
distance function, defined as the distance to the nearest boundary, normally available in practical 
schemes. Applicability of the EM method is limited to model problems with analytical boundary 
shape. 

The more general AM method approximates the EM method by applying the least-square min- 
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imization in a locally constructed coordinate system 


r = f 0 + K?+\v'. (i5) 

The local coordinates, (£', ?/), are constructed using the distance function, which provides infor- 
mation on the closest boundary point. First the distance function is reconstructed at the face center. 
The coordinate vectors at the face center are defined as a unit ^'-directional vector pointing in the 
direction opposite to the closest boundary point and its orthonormal ^'-directional vector. The rf- 
coordinate is taken as the distance from the boundary and the ^'-coordinate is the projection of the 
vector connecting the stencil point with the face center onto the ^'-direction. For the cylindrical 
geometry, the distance function is the shifted radial function, r — R. 

The gradient approximation accuracy for high-F grids of types (I)-(III) from the previous 
study 23 supplemented with the CC-CS and CC-NA-AM results is summarized in Table 5. Con- 
vergence of the maximum gradient errors over all faces is tabulated. Note that large O(Ahg) 
relative errors for the CC-NA scheme occur on high-F grids of type (III) only at the radial faces in 
the gradient component tangential to the face; the errors at other faces and in the gradient compo- 
nent normal to the radial face are small. The large relative errors for the CC-CS scheme occur for 
extremely large T 10 5 at the circumferential and diagonal faces. The least-square gradient itself 
is accurate at these faces; the error is caused by inclusion of the edge gradient in reconstruction. 

Table 5. High-F grids: relative errors of gradient reconstruction 



(I) 

(II) 

(HI) 

NC 

0(h |) 

0(hg) 

0(hg) 

CC-NN 

0(H) 

0(1) 

0(1) 

CC-NN-EM 

o(H) 

0(H) 

0(H) 

CC-NN-AM 

0(h |) 

o(HD 

o(H) 

CC-CS 

o(H) 

O(Ahg) 

O(Ahe) 

CC-NA 

o(H) 

0(H) 

O(Ahg) 

CC-NA-AM 

o(H) 

o(H) 

0(H) 


B. Discretization error convergence 

Discretization errors of cell-centered schemes are compared with the errors of the node-centered 
scheme on refined stretched high- 1 grids of types (III) and (VI). The tests are performed for the 
manufactured solution U = sin(57rr). The computational grids (see Figure 11) are derived from 
background regular cylindrical grids with radial extent of 1 < r < 1.2 and angular extent of 20°. 
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The background grids have four times more nodes in the radial direction than in the circumferential 
direction. The grid-refinement study is performed on grids stretched in the radial direction with a 
fixed maximal aspect ratio A & 1000. The maximal value of parameter T changes approximately 
from 24 to 3. The stretching ratio is changing as [3 — 1.25, 1.11, 1.06, and 1.03. 

* — - NC NC 




(a) Convergence on grids of type (III) (b) Convergence on grids of type (VI) 

Figure 12. Convergence of the L\ -norms of the discretization errors on high-r grids 

Convergence of the Lx -norms of the discretization errors on grids of type (III) is shown in 
Figure 12(a). All tests have been performed stochastically. The plot symbols again indicate the 
mean errors and bars indicate the maximum and minimum errors observed on each scale. As 
expected, error variations observed on grids of the same scale due to grid irregularities are small 
for all schemes and decreasing for smaller scales (larger grids). The errors of the NC, CC-NA, 
CC-CS, CC-NN-AM, and CC-NA- AM solutions converge with second order and are almost over- 
plotted on fine grids indicating the same accuracy per degree of freedom. The errors of the CC-NN 
scheme are significantly higher and converge with first order. 

Convergence of the Li -norms of the discretization errors on grids of type (VI) is shown in 
Figure 12(b). The effective meshsizes of cell-centered and node-centered formulations are much 
closer on mixed grids than on triangular grids. The CC-CS scheme is omitted because its current 
version has been derived only for triangular cells. It is expected that a mixed-element version can 
be derived, but it is not currently available. Note also that the CC-NA scheme may lose stability 
on high-T mixed element grids. On these grids, there are topologies where the node solution is 
averaged from four neighboring cells. The four cell centers involved in such averaging may be 
located on a straight line, thus leading to degeneration. In these (rare) instances, large negative 
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contributions appear on the main diagonals of the full linearization matrix. The scheme may still 
be solved and even provide a reasonable accuracy. The approximately mapped version of the node 
averaging scheme, CC-NA-AM, is always stable. 


IX. Conclusions 

Node-centered (NC) and cell-centered schemes have been compared for finite-volume dis- 
cretization of Poisson’s equation as a model of the viscous flow terms. Among the considered 
schemes, the node-centered scheme has the lowest complexity, i.e., involves the least number of 
neighbors. The cell-centered schemes using least-square face gradient reconstruction have com- 
plexity comparable to that of the node-centered scheme. Complexity of the cell-centered node- 
averaging (CC-NA) scheme is the highest. 

The accuracy comparisons have been made for two classes of tests: the first class is repre- 
sentative of adaptive-grid simulations and involves irregular grids with discontinuous metrics; the 
second class is representative of high-Reynolds number turbulent flow simulations over a curved 
body. All tests have been performed for smooth manufactured solutions. 

For the tests of the first class performed in rectangular geometries on consistently refined irreg- 
ular grids, all schemes demonstrated similar qualities: 

(1) The discretization errors converge with second order and quantitatively are similar for all the 
tested schemes on grids with the same number of degrees of freedom. 

(2) Gradient reconstruction may produce 0(Ah x ) large relative errors on grids of types (IV)- 
(VI); here A is the grid aspect ratio and h x is the larger mesh spacing. 

(3) Truncation errors do not converge, as expected. 

The CC-NA scheme with clipping can fail to approximate gradients and/or to converge to the exact 
solution. Note, however, that the clipping is introduced mainly for stability of the inviscid solution 
and can be eliminated for the viscous terms. 

The tests of the second class have been performed on consistently refined stretched grids gen- 
erated around a curved body, typical of those generated by the method of advancing layers. The 
range of grid parameters has been chosen to enforce significant curvature-induced grid deforma- 
tions, characterized by parameter T. These high-T tests proved to be more discriminating: 

(1) The discretization errors are small and converge with second order for the NC scheme, for 
the cell-centered least-square schemes with a compact stencil (CC-CS) and for approximate 
mapping schemes (CC-NN-AM and CC-NA-AM). The discretization errors for the CC- 
NA scheme also converge with second order, but the scheme may lose stability on mixed 
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grids. The cell-centered nearest-neighbor scheme without mapping (CC-NN) shows first- 
order convergence and the highest level of discretization errors. 

(2) Accurate gradient reconstruction is provided by the NC scheme and the approximate map- 
ping schemes, CC-NN-AM and CC-NA-AM, on all grids. On high-T grids of types (II) and 
(III), the CC-NN scheme generates 0(1) errors in gradient reconstruction. The cell-centered 
CC-NA and CC-CS schemes may produce large relative gradient errors proportional to the 
product of the grid aspect ratio and the larger mesh spacing. 

The major conclusion is that accuracies of the node centered and the best cell-centered schemes 
are comparable at equivalent number of degrees of freedom. 


A. Node-centered discretization on grids of types (II) and (III) in 

rectangular geometries 


In this section, we show that the node-centered (Galerkin) discretization of the Laplacian re- 
duces to a common finite-difference formula for triangular grids for rectangular geometries of 
arbitrary aspect ratio. Consider a set { T t } of triangles/tetrahedra that share a node j. With solution 
values stored at nodes, the Green-Gauss gradient within each cell is given by 


W T 


1 

GQ 


y UiUi, 


(16) 


where D = 2 for triangles, D = 3 for tetrahedra, is the volume of the cell T, {i T } is a set of 
nodes of the cell T, and n, is the scaled inward normal vector of the face opposite to the node % 
whose magnitude is equal to the face area. Then, applying the node-centered discretization as in 
(3) (or equivalently the standard Galerkin discretization), we obtain a discretization at j, 

f AUddj = [ W- nrfT, -- V —O V Vi (n, • nj) , (17) 

■> Te{Tj} T ie{i T } 


where Qj and T y are respectively the dual control volume around j and its boundary, and n J is the 
scaled inward normal vector opposite to the node j in the cell T. For our purpose, it is particularly 
convenient to separate the right-hand side of (17) into two terms: 
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(18) 
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Figure 13. Node-centered stencil on triangular grids. Note that the normals are not scaled. 

The first term contains contributions from the node value, U v and the second term contains contri- 
butions from the neighbors. 


We now restrict our attention to two dimensions and specialize (18) to triangular grids (Figure 


13): 


A UdZli = 



(19) 


where { kj } is a set of neighbors of j and the normals are scaled inward normals defined as in 
Figure 13. This can be written also in terms of angles between faces, 

Uj + \ ( cot + cot Uk ’ ( 2 °) 

k£{kj} 

and is often used to show that the discretization is positive for triangulations with 0l + Or < 7T. 
Consider now a triangular grid shown in Figure 14, which is constructed by inserting diagonals 
into a Cartesian grid. For this particular diagonal splitting, the node 3 does not contribute to the 
discretization (20) because this node is not a neighbor to the node j, and also the nodes 1, 5, 
and 7 do not contribute because the angles Ol and Or are both 90° and therefore the coefficient 
(cot 0 L + cot Or) vanishes. This is, in fact, true for any diagonal splittings', contributions from the 
corner nodes, 1, 3, 5, 7, are always zero either because it is not in the actual stencil or because the 
coefficient vanishes. Observe also that the angles Ol and Or for other nodes are independent of the 
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Figure 14. Triangular grid with arbitrary aspect ratio. 


diagonal splitting, and thus we always have 


cot 9l = cot Or = 


hy 

K 

K 

hr,. 


for nodes 2 and 6, 
for nodes 4 and 8. 


( 21 ) 


Moreover, it is easy to show that the coefficient of Uj is also independent of the splitting. Hence, 
the discretization (20) can be written, for arbitrary splittings, as 
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( 22 ) 

(23) 


This is nothing but a common finite-difference discretization. Therefore, the node-centered (Galerkin) 
discretization of the Laplacian is equivalent to a common finite-difference discretization for tri- 
angular grids for rectangular geometries with arbitrary diagonal splitting and aspect ratio. For 
stretched grids where h x varies in x-direction and/or h y varies in ^/-direction, the comer nodes still 
do not contribute to the discretization. 

In three dimensions, the node-centered discretization does not generally reduce to a common 
finite-difference discretization even for a regular tetrahedral grid which can be constructed from 
a Cartesian grid by dissecting each cube into five or six tetrahedra. It can become equivalent 
precisely to the standard finite-difference discretization in a very special case where a node is 
originally centered at a common 7-point finite-difference stencil. In general, a solution at a node 
will be coupled with almost all solutions at the edge-adjacent neighbor nodes because the dot 
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product (rij • nj) in (18) rarely vanishes. However, actual coefficients for neighbor solutions are 
generally much larger at nodes in a common finite-difference stencil than at other nodes. 
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